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We present a general framework to study stability of the synchronous solution for a hypernetwork 
of coupled dynamical systems. We are able to reduce the dimensionality of the problem by using 
simultaneous block-diagonalization of matrices. We obtain necessary and sufficient conditions for 
stability of the synchronous solution in terms of a set of lower-dimensional problems and test the 
predictions of our low-dimensional analysis through numerical simulations. Under certain conditions, 
this technique may yield a substantial reduction of the dimensionality of the problem. For example, 
for a class of dynamical hypernetworks analyzed in the paper, we discover that arbitrarily large 
networks can be reduced to a collection of subsystems of dimensionality no more than 2. We apply 
our reduction techique to a number of different examples, including a class of undirected unweighted 
hypermotifs of three nodes. 



I. INTRODUCTION 

Much recent work has been devoted to the study of dy- 
namical networks £Q. A few studies have considered the 
dynamics of hypernetworks where the individual units 
are coupled through two or more interaction networks. 
Hypernetworks arise in applications as different as the 
spread of epidemic diseases 0, computer viruses [3], 
game theory [4], social interactions [5], and neural net- 
works formed of both electrical gap junctions and chem- 
ical synapses [B]. In complex adaptive systems, different 
types of couplings usually coexist, including cooperative, 
competitive and symbiotic couplings [7J. 

In this paper, we use the term hypernetwork to in- 
dicate a set of nodes that are coupled through connec- 
tions of different types, with the connections of the same 
type forming a distinct network layer. A similar concept, 
which has been used to describe mainly social systems, is 
that of a multislice or multiplex network [5] , where both 
intra-layer and inter-layer connections are present. Inter- 
dependent networks are usually evoked in the context of 
engineering or technological applications, when the nodes 
in each layer rely on their connections to nodes in other 
layers for their proper functioning (e.g., the coupled func- 
tions of the power grid and computer communication net- 
work studied in [5J). Transportation networks have been 
described as layered networks in [§]. Another definition 
is that of networks of networks, which are evoked in gen- 
eral when connections exist between nodes belonging to 
different networks (see e.g., [TfllITT] ). 

We are interested in the synchronization dynamics 
of hypernetworks. There are many different types of 
synchronization, including complete synchronization |12| - 
[13], phase synchronization [15], lag synchronization [15] , 
group or cluster synchronization |17L 118] , and generalized 
synchronization [19, 20 . For a review of synchronization 
of complex networks, the reader is referred to |21j . Syn- 
chronization of hypernetworks has relavance to the study 
of any system exhibiting multiple types of coupling. For 
example, studying excitation patterns of neural networks 
involves the analysis of both chemical and electrical sig- 



nals between neurons |17) . As another example, studying 
the synchronous motion of entire schools of fish involves 
the analysis of not only the visual cues between fish but 
also the release of chemical signals into the water [2"2" l l23 | . 
In this paper, we focus on complete synchronization of 
hypernetworks. 



The problem of synchronization of dynamical hyper- 
networks has been first studied in (6] , where special con- 
ditions have been considered, for which the problem of 
stability of the synchronous solution can be reduced in 
a low-dimensional form. However, a general framework 
to study stability of the synchronous solution for a dy- 
namical hypernetwork is lacking. In this paper, we will 
present a general approach to obtain a reduction of this 
problem in a low-dimensional form. We will do that by 
looking at lower-dimensional graphs, whose stability will 
characterize that of the original higher-dimensional net- 
work. Moreover, the approach we present can be used to 
find out to what extent the dimensionality of the original 
problem can be reduced. 



Low-dimensional approaches have proved helpful in an- 
alyzing the dynamics of networks of coupled dynami- 
cal systems. Examples include (i) the stability of the 
synchronous evolution for networks of coupled oscillators 
[II mini [HHg, (ii) the stability of the consensus state 
in networks of coupled integrators , (iii) the stability 
of discrete state models of genetic control [37], and (iv) 
the stability of strategies in networks of coupled agents 
playing a version of the prisoner's dilemma |28j . In this 
paper, we are interested in studying the stability of the 
synchronous solution for a dynamical hypernetwork and 
we show that reducing the stability problem in a lower- 
dimensional form is an available approach. 
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II. MODEL 

We consider a dynamical hypernetwork described by 
the following system of coupled differential equations: 

M N 

xi(t) = F( Xi {t)) + Y,Y, A S )H ^(t - **)). (!) 

fe=l j=l 

i = 1, ...N, where Xi(t) is the m-dimensional state of sys- 
tem i at time t, i = 1, N, the function F : R m -> i? m 
determines the dynamics of each individual system when 
uncoupled; H k : i? m — > R m are arbitrary coupling func- 
tions, k = 1,...,M, Tfc > is the time-delay associated 
with the coupling function H k , k = 1, ...,M. The entries 
of the matrix A^ = {A, ( *°} are such that ^ if 
node j is coupled to node i through the coupling func- 
tion H k and Ay = otherwise. Moreover, we require 
that J2j Ajty — a k , i.e., the sum of the entries along the 
rows of each matrix Aw is constant and is equal to a k 
[55] , Under this assumption, system ([T]) allows the fol- 
lowing synchronous solution, 

xi(t) =x 2 {t) = ... = x N {t) =x s (t), (2) 

obeying 

x. (t) = F(x s (t))+J2 akH k {x s [t-r k )) = F(x s (t) ). (3) 
fc 

Now, by replacing the function F with the function F in 
([3]), we obtain, 

M N 

x l (t) = F(x l (t)) + Y / Y, L ^ Hk ^ t - T ^> ( 4 ) 
fc=U=i 

where each matrix L( fc ' = {XV/}, with zi*) = A,-!/ — 

(5,.,a fc has the property that Y^j ^ij = 0, k = 1, ...,M, 
i.e., the sum of the elements in each row is zero, and 
following a common convention [T], we refer to such ma- 
trices as Laplacian matrices. 

In order to study stability of the synchronous solution, 
we linearize Eqs. Q about ([2|, obtaining, 

S Xi (t) =DF(x s {t))5 Xi (t) 

M N ,r\ 

+ 4?DH k (x s (t - T k ))S Xj (t - Tfc), 

k=l 3=1 

or equivalently, in matrix form, 

SX(t) =I N ® DF{x s (t))8X[t) 

M (6) 
+ J2 L (k) <8 DH k (x s (t - r k ))5X{t - r k ) 

k=l 

where the mTV-dimensional vector X(t) = 
[xi(t) T , X2(t) T , XN(t) T ] T and we have used the 



symbol (g> to indicate the Kronecker product or direct 
product. 

For the sake of simplicity, in what follows, we focus on 
the case that M = 2, for which ^ can be rewritten, 

SX(t) =I N ® DF(x s (t))5X(t) 

+L (1) QDH^Xsit-TiftSXit-n) (7) 
+L (2) ® DH 2 (x s (t - t 2 ))5X(< - r 2 ). 

Note that ([7]) is an mTV-dimensional system, in the 
sense that it is described by mN coupled state variables. 
[57] . Our goal will be to reduce the dimensionality of the 
system (JtJ) , by decoupling the stability problem Q into a 
set of lower-dimensional problems, each one independent 
of the others. 

It was shown in Ref. [5] that there are three cases for 

which the mN dimensional problem Q can be reduced 

to a set of (N — 1) 2m-dimensional problems. These 

three cases are: (i) the Laplacian matrices and 

commute; (ii) one of the two networks, say k — 2, is 

(2) 

unweighted and fully connected, i.e., L\- — c for i 7^ 

j, = -c(N - 1), i = 1,...,N; (iii) one of the two 
networks, say k = 2, is such that the coupling strength 
from node i to node j is a function of j but not of i, i.e., 

Lf> r, for i •• ./. ii 2) = - 1- •••• N. _ 

However, if none of the three above conditions is sat- 
isfied (and each of the conditions (i),(ii),(iii) generally 
do not occur if the coupling strengths of the networks 
are arbitrarily chosen), such a reduction is not possi- 
ble. In what follows, we will extend the results of [5] 
with the goal of reducing the original stability problem 
to a set of n subproblems of maximum dimension a, with 
2 < a < N, depending on the properties of the matrices 
I», k = 1,...,M. 

In general, we will be interested in addressing the fol- 
lowing algebraic problem: given the set of ./V-square real 
matrices L = {L^jL' 2 ',...,^'}, find the finest simul- 
taneous block-diagonalization (SBD) of C. The prob- 
lem consists in finding an invertible matrix P, such that 
p-ijjfyp = ©« =1 Bj ) where the symbol © denotes the 

direct sum of matrices, Bj , Bj , 2?™ are square matrix 
blocks of dimension bj, and J2^ =1 bj — N. The diago- 
nalization is said to be the finest if the maximum block 
dimension b max = max™ =1 bj is minimal with respect to 
the choice of P. The problem can either be solved ex- 
actly or in a parameterized (approximate) form, where 
for a given parameter e > 0, P^L^P = ®j =1 Bj 

i = 1, ...,M, and the N dimensional matrices E^ % > , 
i = 1, M, are of order e in magnitude. 

Hereafter, we briefly review an algorithm for SBD of 
sets of matrices. Such a block-diagonal decomposition is 
not unique in general and naturally we are interested in 
finding the matrix P that provides the finest decomposi- 
tion. 

Instead of trying to tackle the problem directly, the 
approach in [H],[3n] aims at finding a basis that di- 
agonalizes the ^-algebra associated with the algebra 



3 



generated by C. This corresponds to finding a ma- 
trix U that simultaneously commutes with the matrices 
LW, L( 2 \ L,( M \ i.e., that simultaneously satisfies the 
following sets of equations: 



UL (i) _ L {i)jj = Q 



(8) 



i = l,...,M. 

The steps of the algorithm are described in what fol- 
lows: 

(i) Let O® be the iV 2 -matrix =I N ® L w - L (l) ® 
In- 

(ii) Construct the matrix S = {i)T . 

(iii) Let y be any AT 2 -vector in the null subspace of 
the matrix S. The iV 2 -vector u can be subdivided in N 
vectors of dimension N as follows, u — [iifjU^, ..■,u 1 N ] T . 

(iv) Obtain U as the matrix whose columns are u\, 

U 2 ,..., U N . 

(v) Finally, P can be constructed as the matrix whose 
columns are the eigenvectors of U. 

In certain situations (when for example, a satisfac- 
tory SBD reduction is not available), we might be in- 
terested in finding a parameterized (approximate) SBD, 
which can be formulated as follows. Given a parame- 
ter e > 0, find the invcrtible iV-square matrix P such 
that P- X L^P = ffiJUBj i = 1,...,M, and the 

N dimensional matrices E^ l \ i = 1,...,M, are of order 
e in magnitude. An error controlled version of the SBD 
algorithm can be found in Ref. (16). 

Now consider problem Q . Suppose we have been able 
to find the finest SBD for C = {L^\ L^>} and that this 
is provided by the invertible matrix P. Then we left- 
multiply both sides of Eq. (|7| by P^ 1 ® I m and by using 
the change of variables rj(t) = P^ 1 ® I m SX(t), we can 
rewrite Q as follows, 

j){t) = I N ® DF(x s (t))ri{t) 

+ ( ©7=i B [ p) <g> DH 1 {x s {t ~ n))t](t - n) (9) 
+ ( ©7=i Bf ] ) ® DH 2 (x s (t - r 2 )) V (t - r 2 ). 

It is easy to see then that the system of equations ([9| 
can be decomposed into the following n subsystems, 

rjjit) =I hj ® DF(x s {t))r jj (t)+ 



D 



(2) 



DH 2 (x s (t-T 2 ))r) j {t-T 2 ), 



j = 1, ...,n, where each vector rjj(t) has dimension mbj 
and evolves independently from the others. Moreover, 
J2j bj — N. Each of these subsystems is forced by the 
synchronous evolution x s (t), which obeys Eq. |3|. Thus 
the original miV-dimcnsional problem has been reduced 
to n lower-dimensional problems, each with dimension 
m(bi + l),m(b 2 + 1), ...,m(b n + 1). Also, our proposed 
goal has been achieved with a — (b max + 1), where 2 < 
a < N. Moreover, this reduction is the finest, in the 



sense that it is not possible to obtain another reduction 
in an (a— l)-dimensional form, provided that the original 
block-diagonalization was the finest. 

We note that by construction, the matrices L^ k \ k = 
1, M, have a zero eigenvalue with associated eigenvec- 
tor [1, 1, 1] T . Hence, when obtaining the finest SBD, 
there must be a 1-dimensional subsystem, indexed j = 1, 



for which B 



(AO 



= 0, k = 1, M, yielding, 
77! (t) = DF(i,(t))»n(t). 



(11) 



We note that this one equation is associated with per- 
turbations lying in the direction parallel to the synchro- 
nization manifold (given by the eigenvector [1, 1, 1] T ) 
and therefore it is irrelevant in determining transversal 
stability of the synchronous solution. 

For a generic subsystem of dimension D, we may want 
to identify a minimal set of parameters (pi,p 2 , ■■■iPm) 
that characterize stability. In general, it can be shown 
that the minimum number of parameters is td — D 2 + 1 
(see Sec. HA). However, for some specific cases, it is 
possible to parameterize a D-dimensional subsystem by 
using less than Tn parameters. 



A. The particular case that 6„ 



B) ' ® DHx{x„{t - Tx))rn(t - n)+ (10) p| becomes 



We now look at Eq. 10 and consider the particular 
case that b max = 2, i.e., a — 3. We further assume that 
the simultaneous block-diagonalization yields /i blocks of 
dimension 1 and v blocks of dimension 2, with (fi + 2u) = 
(N — 1). The problem that we want to address in this 
section is described below. 

Consider all the pairs of matrices C = {L^,L^} 
such that a simultaneous block-diagonalization can be 
achieved with b max = 2. For blocks of dimension D = 
{1, 2}, we aim at finding a reduction of the stability prob- 
lem in a parametric form in the scalar parameters 
(pi,P2, ■■■,Pr D ), such that 7\d is minimal. In what fol- 
lows, we independently address this problem for blocks 
of dimension D — 1 and blocks of dimension D = 2 and 
we show that r\ = 2 and r 2 = 5. We also obtain a general 
relation between rjj and the block dimension D. 

It is easy to see that for blocks of dimension 1, Eq. 



fij(t) =DF(x s {t)) Vo {t) 

+Bf ) DH 1 {x s {t-T 1 ))i lj {t-T 1 ) 
+Bf ) DH 2 {x s {t-T 2 ))r l3 {t-T 2 ) 



(12) 



(i) 



j = l,...,/i, where rjj{i) has dimension m and B 

/ 2 \ 

and B)j are two scalar (eventually, complex) parame- 



ters. Each subsystem (12 1 is parametrized by the pair 



{Bf\Bf ] ). Hence, r x = 2 
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For blocks of dimension 2, Eq. (10 1 becomes, 

r)j(t) =I 2 ®DF(x s (t))r] j (t) 

DHiixsit-nVrjjit-n) (13) 



-D 



(i) 



+Bf ® DH 2 (x s (t - t 2 )) % -(< - r 2 ) 



II M 




j = 1. ^, where has dimension 2m and sj 1 ' and 

(2) 

Bj are two square matrices of dimension 2. We note 
that it is possible to further diagonalize either one of the 



two matrices B^p or B^'; without loss of generality, we 



? (2). 




diagonalize Bf\ obtaining Bf ] = WjKjW' 1 . By pre- '{ FCgQ^QfCG ) 



multiplying each block (13) by W, 1 ® I m , we obtain, 



Q(t) =I 2 ®DF(x s {t))(j{t) 

+Aj ® DHi(x s (t - n))Cj(t - n) 



(14) 



+Qj <g> DH 2 {x s (t - T 2 ))C,j{t - r 2 ), 

ImVjit), Aj is the 



j = where (*) = Wj 1 

following diagonal matrix, 



A; 







X 2 



and the matrix Qj = Wj 1 B^Wj is the following, 



qf qf 
S 1 q¥ 



(15) 



(16) 



We see that for each block j = l,...,u, stabil- 
ity depends on the following set of sca lar parameters 
(Aj, Af , qf, qf,qf > if )• From < |13|15|16[> we see that Eq. 
( 14 ) can be decomposed into the following two coupled 



equations, 

C)(t) =DF(x s (t)){}(t) + x)DH 1 (x s (t - n ))cj(t - n) 

+qfDH 2 (x s (t-T 2 ))(}(t-T 2 ) 
+qfDH 2 (x s (t-T 2 ))(](t-T 2 ), 

C?(t) =DF(x s (t)){ 2 (t) + xpH^t - n ))$(t - n ) 



+q 21 DH 2 (x s (t-r 2 ))(}(t 



T2 



+q 22 DH 2 (x s (t-T 2 ))C 2 (t-T 2 ), 



where the vector C, 



•3 ~ i>j 

substitution, qf( 2 (t) — >• Cf (^)i we can rewrite (17), 



Kf(t)Xf(t)] T . 



(17) 

Now, with the 



#(t) =DF(x a {t))$(t) + X)DH 1 (x s (t - n))Cj(t - ti) 

+ (Z f J Di/ 2 ( 2 ; s (t-r 2 ))C J 1 (i-r2) 

+ J Di? 2 (x s (i-T 2 ))C J 2 (t-r 2 ), 
C 2 (i) =DF(x s (t))g(t) + X 2 DH 1 (x s (t - n))C 7 2 (t - n) 



+qfqfDH 2 {x s {t - r 2 ))Cj(t - r 2 ) 
+gfDfl-2(x 4 (t-75))C?(t-75). 



(18) 



CP 

CP 
CP- 



-N-3 



FIG. 1: [Color Online] A special class of hypernetwork con- 
figurations. These hypernetworks shown on the left can be 
neatly reduced into the (N — 2) subsystems shown on the 
right. The reduction yields a single two-node system, followed 
by (N — 3) one-node systems. Stability of the hypernetworks 
on the left corresponds to that of the lower-dimensional sub- 
systems on the right. 



Thus each subsystem (18), j = l,...,v, is described 
by the following set of r 2 = 5 scalar parameters 
(Aj, Xp qf ,qfqf ,q 22 ). It follows that for 2-dimensional 
subsystems, r 2 = 5. 

We conclude that each one-dimensional subsystem j — 
l,...,/i can be associated with a master stability func- 
tion M 1 (Bj 1 \ Bj ) which returns the maximum Lya- 
punov exponent of Eq. (12) as a function of the pair 



(B^ ,B)j ). Also, each two-dimensional subsystem j — 
l,...,u can be associated with a master stability func- 
tion 7Vf 2 (A], X 2 , qf, qfqf, q 22 ) which returns the max- 



imum Lyapunov exponent of Eq. (17) as a function of 
the 5-tuple (A], A^, qf , qf qf , q 22 ) . Once the master sta- 
bility functions A4 1 , M 2 are known, stability of the syn- 
chronous solution for a generic dynamical hypernetwork 
(described by Eq. (1) with M = 2) that allows a simulta- 
neous block-diagonalization with b max = 2, can be deter- 
mined by knowledge of the pairs (Bj , Bf 1 ) for blocks 
of dimension 1 and of the 5-tuples (Aj, A|, qf , qf qf , q 22 ) 
for blocks of dimension 2. 

By extending the above reasoning to subsystems of 
higher dimension D, it can be shown that rrj = D + 1. 



III. NUMERICAL EXAMPLES 

The left hand side of Fig. [T] shows a special class of hy- 
pernetworks for which the stability problem can be con- 
veniently reduced by using SBD. The class contains all 
hypernetworks made from two identical fully connected 
graphs (FCG), each of size ^, that are connected to one 
another only by a single alternative connection. The pa- 
rameter a is the coupling strength of all the connections 
inside each FCG and the parameter b is the coupling 
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strengths of the alternative connection. The associated 
Laplacian matrices do not commute unless either a = 
or 6 = 0. As can be seen on the right hand side of Fig. 
[T] we discover that a hypernetwork in this configuration 
can always be reduced to a collection of subsystems of 
dimensionality no more than 2, that is, one subsystem 
of dimension 2 and (N — 3) identical subsystems of di- 
mension 1 (we neglect the one subsystem that is associ- 
ated with perturbations parallel to the synchronization 
manifold). The reduction in this form becomes partic- 
ularly advantageous when the dimension of the original 
hypernetwork iV is large. From the SBD decomposition 
(described in Sec. II), we obtain that the parameter d in 
the figure depends upon the size of the hypernetwork, i.c, 
d=a(f -1). 

As an example, we consider the hypernetwork shown 
on the left hand side of Fig. [l]of dimension N = 6. Its 
dynamics is described by the set of Eqs. where each 
individual node obeys the equation of the Lorenz chaotic 
system, for which m = 3, x(t) = (xi(t), X2(t), x-z(t)) T , 



F(x) = 



W[x 2 (t) - Xl {t)] 
xi(t)[28-aj 3 (t)] -x 2 {t) 
xi{t)x2{t)-2x s (t) 



(19) 



fli(a:(i)) = [0,x 2 (t),0] T , H 2 (x(t)) = [ Xl (t),0,x 3 (t)} T , 
and t\ = T2 = 0. The adjacency matrices and A^ 
correspond respectively to the black [black] and gray [red] 
connections of the N — 6 hypernetwork in Fig. [Tjand are 
defined as follows: A^) = Af) = a if (i - 0) x (j -6)>0, 
otherwise, where 9 = A\j — Ajf = except for 

, j = j*), with 1 < i* < ^ and 

" (fc) -(4 fc) -%E,4 fc) ).fc 



the one pair (i = » ,j — j j, «™i a „ ^ ., 

From the SBD procedure, we obtain that stability of 
the original 7Vm-dimensional system can be reduced to 
that of two lower-dimensional systems, one of dimension 
m and one of dimension 2m (see Fig. [I]). We indepently 
compute the maximum Lyapunov exponent (MLE) for 
both these systems. For the m-dimensional subsystem, 
we find that the condition for stability is that a^- > 2.29. 
For the 2m-dimensional subsystem, we record the maxi- 
mum Lyapunov exponent as a function of the pair (a, b) , 
for the specific case of N = 6. The results of our numer- 
ical computations are summarized in the upper plot of 
Fig. [2j where the light gray [yellow] area corresponds to 
the region of the (a, b) plane for which the MLE of the 
m-dimensional subsystem is negative and the dark gray 
[gray] area corresponds to the region of (a, b) plane for 
which the MLE of the 2m-dimensional system is nega- 
tive. Note that for this case, the intersection coincides 
with the gray area. 

In order to test our low-dimensional predictions, we 
numerically integrate Eqs. ([!]) from an initial condition 
close to the synchronization manifold. For each run, we 
monitor the average synchronization error E, 




E(t) = (NAt) 



N 



t+At 



\\xi(r) - x(T)\\dr, (20) 



FIG. 2: [Color Online] The upper plot shows the areas of the 
(a, b) plane corresponding to a negative MLE for both the Tri- 
dimensional subsystem (colored in light gray [yellow]) and the 
2m-dimensional subsystems (colored in dark gray [gray]). We 
expect stability of the original hypernetwork in the intersec- 
tion of the light gray [yellow] and dark gray [gray] areas. The 
lower plot shows simulations of the full high-dimensional hy- 
pernetwork, N = 6; we plot in gray the area of the (a, fe)-plane 
for which the synchronization error E(t) decreases steadily 
below 1% of its initial value. 



where x(t) = N^ 1 Y^iLi x i(t) an d ll£ll indicates the Eu- 
clidean norm of the vector £. The lower plot of Fig. [2] 
shows the area of the (a, 6)-plane for which E(t) is ob- 
served to decrease steadily below 1% of its initial value. 
As can be seen from the upper and lower plots of Fig. 
[2] there is very good agreement between the dynamics of 
the original hypernetwork and its low-dimensional coun- 
terpart. 

We also considered the case shown in Fig(3] that the 
two FCG graphs are connected by R alternative connec- 
tions rather than 1, each one with associated strength b, 
and such that the endpoints of these R connections never 
coincide in the same node. How is the stability of this hy- 
pernetwork going to be characterized? In what follows, 
we restrict our attention to the case that 1 < R < N/2. 
By applying the SBD procedure, we discover that, as can 
be seen from Fig. [3j a hypernetwork in this configuration 
can always be reduced to a collection of subsystems of di- 
mensionality no more than 2, that is, one subsystem of 
dimension 2 and (N — 3) subsystems of dimension 1 (we 
neglect the one subsystem that is associated with pertur- 
bations parallel to the synchronization manifold). The 



2b] 



aN 
~~2~ 



FIG. 3: [Color Online] On the left: a hypernetwork formed 
of two identical fully connected graphs, each of size $ , con- 
nected to one another by a set of R alternative connections, 
1 < R < N/2. The endpoints of the alternative connections 
never coincide in the same node. The case that R — 1 corre- 
sponds to that studied above in Fig. [T] These hypernetworks 
shown on the left can be neatly reduced into the (N — 2) 
subsystems shown on the right. The reduction yields a single 
two-node system, followed by (N—3) one-node systems, which 
are divided into two distinct groups of (R— 1) and (N — 2 — 7?) 
identical subsystems. Stability of the hypernetworks on the 
left corresponds to that of the lower-dimensional subsystems 
on the right. 



one subsystem of dimension 2 depends on the number of 
alternative connections R, with the parameters shown in 
the figure d! = — R) and a' = aR. Moreover, for 
this more general case, as can be seen on the right hand 
side of Fig. |3j the remaining (N — 3) subsystems of di- 
mension 1 are divided in two distinct groups of (R — 1) 
and (N — 2 — R) identical subsystems. Then, in order 
to characterize stability of hypernetworks in this config- 
uration, we have to consider stability of all the three dif- 
ferent types of subsystems that arise from the reduction. 
Stability would occur in a region in the parameter space 
which is the intersection of the three resulting stability 
regions. The reduction would still be very significant for 
large enough values of N. 

As another example, we consider the hypernetwork 
shown on the left hand side of Fig. [4] From the figure we 
see that the hypernetwork is composed of M = 2 different 
networks, each one associated with a different coupling 
function Hk (to be defined in what follows). Hence, the 
linearized problem can be cast exactly in the form of Eq. 
(7) with M = 2, where the two Laplacian matrices are 
as follows, 




(21) 



FIG. 4: [Color Online] On the left, a hypernetwork formed 
of two interaction graphs. Each node represents a linearized 
system. On the right, we have obtained a reduction into two 
lower dimensional systems: one of b\ — 2 nodes and another 
one formed of 62 = 1 node. 



associated with the black [black] connections in the figure 
and 



—a 


a 





a 


—a 











-b 








b 



(22) 



associated with the gray [red] connections in the figure. 

The two matrices do not commute unless b = a. There- 
fore, we consider the case b ^= a. The dynamical hyper- 
network is described by the set of Eqs. 0, where each 
individual node obeys the equation of the Lorenz chaotic 
system, for which m — 3, x(t) — (xi(t) , x 2 (t) , x 3 (t)) T , 



F(x) = 



10[x 2 (t) - x^t)} 
xt(t)[28 - x 3 (t)} -x 2 {t) 
xi{t)x 2 {t) - |aj s (*) 



(23) 



Ht(x{t)) = Mi),0,0] T , H 2 (x(t)) = [0,0,x 3 {t)} T , and 
n = r 2 = 0. 

By using the procedure described in [29j . we find a 
matrix P, 



P = 



-0.5000 -0.5000 -0.4330 -0.5590 

-0.5000 -0.5000 0.4330 0.5590 

-0.5000 0.5000 -0.5590 0.4330 

-0.5000 0.5000 0.5590 -0.4330 



(24) 



that simultaneously block-diagonalizes L^> and L^ 2 \ By 
using P, we obtain that the original 4m-dimensional sys- 
tem can be decomposed in two m-dimensional systems in 
the blocks (0, 0) and (—2a, 0) and in one 2m-dimensional 
system. The subsystem associated with the pair (0, 0) 
corresponds to perturbations parallel to the synchroniza- 
tion manifold and as such is irrelevant in determining 
transversal stability of the synchronous solution. 
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a negative MLE for the m-dimensional subsystem is col- 
ored in light gray [yellow] and the area associated with a 
negative MLE for the 2m-dimensional subsystem is col- 
ored in dark gray [gray] . We expect stability of the orig- 
inal hypernetwork in the intersection of the light gray 
[yellow] and dark gray [gray] areas in the f. 

In order to test our predictions, we numerically in- 
tegrate from an initial condition close to the synchro- 
nization manifold the equations of the dynamical hy- 
pernetwork (1), with M — 2, the function F given in 
@, H^xit)) = [s 1 (t) J J 0] r > H 2 (x(t)) = [0,Q,x 3 (t)] T , 
Ti = T2 = 0, and the two Laplacian matrices and 
£,( 2 ) given in Eqs. (21 1 and (22). For each run, we mon- 
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itor the average synchronization error E, defined in Eq. 
(20 1. The lower plot of Fig. [5] shows the final synchro- 
nization error E versus b for a = 4, which converges 
to zero in the stability interval predicted by the lower- 
dimensional analysis. 



FIG. 5: [Color Online] The upper plot shows the areas of 
the (a, b) plane to which corresponds a negative maximum 
Lyapunov exponent (MLE) for both the m-dimensional sub- 
system on the right-hand side of Fig. 2] (colored in light gray 
[yellow]) and the 2m-dimensional subsystems on the right- 
hand side of Fig. [4] (colored in dark gray [gray]). We expect 
stability of the original hypernetwork in the intersection of 
the light gray [yellow] and dark gray [gray] areas. The lower 
plot shows the synchronization error E versus b for a = 4, 
which converges to zero in the stability interval predicted by 
the lower-dimensional analysis. 



By further diagonalizing the one 2m-dimensional sub- 
system, we obtain that this can be recast into the form, 



Cj(t) =l2^DF(x s (t))Cj(t) 

+Aj ® DH 1 (x s (t - ri))0(f - n) 
+Q 3 <g> DH 2 {x s (t - r 2 ))0(f - t 2 ) 



(25) 



with 



and 



Ai = 





-2a 



Qj = 



-(a + b) 1 
(b -a) 2 -(a + b) 



(26) 



(27) 



The procedure to obtain Eq. (25) is illustrated in Sec. 
IIA, compare with Eq. ( 14 ) therein. The right hand side 
of Fig. [4] shows the two lower-dimensional subsystems in 
which the stability problem has been reduced. We ob- 
serve that for this specific problem, stability of both the 
m and 2m-dimensional subsystems can be conveniently 
parameterized in the pair (a, 6). The upper plot of Fig. 
[5] shows the sign of the maximum Lyapunov exponent 
(MLE) associated with both the m and 2m-dimensional 
subsystems in the (a, b) plane. The area associated with 



A. Dynamical Hypermotifs 

As a further application of our theory, we have consid- 
ered the synchronization of dynamical hypermotifs. Mo- 
tifs were introduced in Ref. |31j as recurrent patterns 
of interconnections occurring in complex networks. Syn- 
chronization of small network motifs has been studied 
in [32 Here, we are interested in hypermotifs, i.e. 
motifs with multiple types of coupling. 

In particular, we have considered the class of all the 
possible N = 3-node unweighted undirected hypermotifs 
with M = 2 connection types. We have assumed all the 
connections to have unitary weights. We have obtained 
a list of 6 different such hypermotifs, excluding those for 
which = and those obtained from one of the 
motifs in the list by interchanging the matrix A^' with 
the matrix A^. 

Figure 6 shows all of these 6 hypermotifs, labeled as 
A-F. We have found that for the hypermotifs D-F the 
Laplacian matrices associated with A^ and A^ com- 
mute. Instead, for the hypermotifs A-C, we have ap- 
plied the SBD procedure to reduce them in their lower- 
dimensional form. Their lower-dimensional counterparts 
are also shown in Fig. 6 (center column). The study of 
more elaborated hypermotifs (e.g., with direct connec- 
tions or with more than 3 nodes) is beyond the scope of 
this paper. 



IV. CONCLUSIONS 

In this paper, we introduced a general framework to 
study stability of the synchronous solution of a dynam- 
ical hypernetwork by means of a dimensionality reduc- 
tion strategy. For any set of arbitrarily chosen coupling 
matrices, we are able to obtain the finest SBD (simulta- 
neous block diagonalization) and to evaluate stability of 
the synchronous solution based on that. Under certain 



aQ— O 

b 

bCRO 




\ ♦** 0.75 



■0.5 



-2i 




-1.5 



\ -2!?? 

0.75 







FIG. 6: [Color Online] A-F are all the possible unweighted undirected hypermotifs with N=3 nodes and M=2 connection types, 
excluding those for which A*- 1 ' = A' 2 ' and those obtained from A-F by interchanging the matrix A' 1 ' with the matrix A^ 2 K 
All the connections in A-F have associated unitary weights. For the hypermotifs A-C we use the SBD procedure to reduce 
them in their lower-dimensional form (the lower dimensional graphs are those in the center column). For hypermotifs D-F, the 
Laplacian matrices corresponding to A^ 1 ' and A*- 2 ' commute. Hence, their dynamical reduction is not shown. 



conditions, this technique may yield a substantial reduc- 
tion of the dimensionality of the problem. For example, 
for a class of dynamical hypernetworks analyzed in this 
paper, we discovered that arbitrarily large networks can 
be reduced to a collection of subsystems of dimensional- 
ity no more than 2. Other times the reduction may be 
less significant. 

We have applied our reduction techique to a num- 



ber of different examples, including small undirected un- 
weighted hypermotifs of 3 nodes. An important advan- 
tage of the SBD decomposition is that it can be used to 
find out to what extent the dimensionality of the original 
problem can be reduced. The study of synchronization 
of large arbitrary dynamical hypernetworks is the subject 
of ongoing investigations. 



I 

The authors are indebted to Jens Lorenz for insightful discussions. 
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Note that the constant-row-sum condition (i.e., that the 
sum of the rows of the matrix A(k) is constant and equal 
to a k ) is more general than the zero-row-sum condition, 
usually considered for complete synchronization [12H14] . 
Even when the constant-row-sum condition is not met, 
its satisfaction can be dynamically obtained by means of 
an adaptive strategy 34, 35 . 

However, if either ti > or T2 > the number of initial 
conditions that are needed to describe the evolution of 
the system can be much larger as for each i = 1, ...,N, 
knowledge of Xi(t) over a time interval of length T max = 
maxr, is required. 



